% The profit function 

function y=pij(hj,hje,heqm,param)
    
J=param(1); c=param(2);
Su=param(3); 

y=(-1).*c.*hj+(2+heqm).^2.*hj.*(1+hj).^(-1).*(1+hje).^2.*(4+2.* ...
  heqm.*J+hje.*(4+heqm+heqm.*J)).^(-1).*(4.*heqm.*(1+heqm).*((-1)+J) ...
  +hje.*(4+(-3).*heqm.^2+4.*heqm.*(1+heqm).*J)+hje.^2.*(4+heqm.*(3+ ...
  J+heqm.*J))).^(-1/2).*Su;

end